#define KOVASZNAY_LAMBDA(NU)                                    \
    ( 1./(2*NU) - vf::sqrt( (1./(4*NU*NU) ) + 4*M_PI*M_PI ) )   \
    /**/
///////////////////////////////////////////////////////////////////////////////////////
#define KOVASZNAY_VELOCITY(NU)                                          \
    vf::vec( 1.0 - vf::exp( KOVASZNAY_LAMBDA(NU)*vf::Px() )*vf::cos(2*M_PI*vf::Py()), \
             (KOVASZNAY_LAMBDA(NU)/(2*M_PI))*vf::exp(KOVASZNAY_LAMBDA(NU)*vf::Px() )*vf::sin(2*M_PI*vf::Py()) ) \
    /**/
#define KOVASZNAY_PRESSURE(NU)                                          \
    (1.0-vf::exp(2*KOVASZNAY_LAMBDA(NU)*vf::Px()))/2.0                  \
    /**/
///////////////////////////////////////////////////////////////////////////////////////
#define KOVASZNAY_DU1DX(NU)                                             \
    (-KOVASZNAY_LAMBDA(NU)*vf::exp( KOVASZNAY_LAMBDA(NU)*vf::Px() )*vf::cos(2*M_PI*vf::Py() ) ) \
    /**/
#define KOVASZNAY_DU1DY(NU)                                             \
    (2*M_PI*vf::exp( KOVASZNAY_LAMBDA(NU)*vf::Px() )*vf::sin(2*M_PI*vf::Py() ) ) \
    /**/
#define KOVASZNAY_DU2DX(NU)                                             \
    (vf::pow(KOVASZNAY_LAMBDA(NU),2)/(2*M_PI))*vf::exp( KOVASZNAY_LAMBDA(NU)*vf::Px() )*vf::sin(2*M_PI*vf::Py() ) \
    /**/
#define KOVASZNAY_DU2DY(NU)                                             \
    (KOVASZNAY_LAMBDA(NU)*vf::exp( KOVASZNAY_LAMBDA(NU)*vf::Px() )*vf::cos(2*M_PI*vf::Py()) ) \
    /**/
#define KOVASZNAY_D2U1DXX(NU)                                           \
    (-vf::pow(KOVASZNAY_LAMBDA(NU),2)*vf::exp(KOVASZNAY_LAMBDA(NU)*vf::Px())*vf::cos(2*M_PI*vf::Py())) \
    /**/
#define KOVASZNAY_D2U1DYY(NU)                                          \
    4*M_PI*M_PI*vf::exp(KOVASZNAY_LAMBDA(NU)*vf::Px())*vf::cos(2*M_PI*vf::Py()) \
    /**/
#define KOVASZNAY_D2U2DXX(NU)                                           \
    (vf::pow(KOVASZNAY_LAMBDA(NU),3)/(2*M_PI))*vf::exp(KOVASZNAY_LAMBDA(NU)*vf::Px())*vf::sin(2*M_PI*vf::Py()) \
    /**/
#define KOVASZNAY_D2U2DYY(NU)                                           \
    (-2*M_PI*KOVASZNAY_LAMBDA(NU)*vf::exp(KOVASZNAY_LAMBDA(NU)*vf::Px())*vf::sin(2*M_PI*vf::Py()) ) \
    /**/
#define KOVASZNAY_D2U1DXY(NU)                                           \
    2*M_PI*KOVASZNAY_LAMBDA(NU)*vf::exp(KOVASZNAY_LAMBDA(NU)*vf::Px())*vf::sin(2*M_PI*vf::Py()) \
    /**/
#define KOVASZNAY_D2U2DXY(NU)                                           \
    vf::pow(KOVASZNAY_LAMBDA(NU),2)*vf::exp(KOVASZNAY_LAMBDA(NU)*vf::Px())*vf::cos(2*M_PI*vf::Py()) \
    /**/
///////////////////////////////////////////////////////////////////////////////////////
#define KOVASZNAY_GRAD_VELOCITY(NU)                                     \
    vf::mat<2,2>(KOVASZNAY_DU1DX(NU),KOVASZNAY_DU1DY(NU),               \
                 KOVASZNAY_DU2DX(NU),KOVASZNAY_DU2DY(NU) )              \
    /**/
    /**/
#define KOVASZNAY_GRAD_P(NU)                                            \
    vf::vec(-KOVASZNAY_LAMBDA(NU)*vf::exp(2*KOVASZNAY_LAMBDA(NU)*vf::Px()), \
            vf::cst(0.) )                                               \
    /**/
///////////////////////////////////////////////////////////////////////////////////////
#define KOVASZNAY_SIGMA_VELOCITY(NU)                                    \
    vf::mat<2,2>(2*KOVASZNAY_DU1DX(NU), KOVASZNAY_DU1DY(NU)+KOVASZNAY_DU2DX(NU), \
                 KOVASZNAY_DU2DX(NU)+KOVASZNAY_DU1DY(NU), 2*KOVASZNAY_DU2DY(NU) ) \
    /**/
#define KOVASZNAY_SIGMA_PRESSURE(NU)                \
    (-KOVASZNAY_PRESSURE(NU)*vf::mat<2,2>(vf::cst(1.),vf::cst(0.),vf::cst(0.),vf::cst(1.))) \
    /**/
#define KOVASZNAY_SIGMA(NU)                                     \
    (KOVASZNAY_SIGMA_PRESSURE(NU) + NU*KOVASZNAY_SIGMA_VELOCITY(NU))    \
    /**/
///////////////////////////////////////////////////////////////////////////////////////
#define KOVASZNAY_DIV_SIGMA(NU)                                         \
    (-KOVASZNAY_GRAD_P(NU) + NU*vec( 2*KOVASZNAY_D2U1DXX(NU) + KOVASZNAY_D2U1DYY(NU) + KOVASZNAY_D2U2DXY(NU), \
                                     KOVASZNAY_D2U2DXX(NU) + KOVASZNAY_D2U1DXY(NU) + 2*KOVASZNAY_D2U2DYY(NU) ) ) \
    /**/

#if 0
#define FORCE_VOLUMIQUE                                                 \
    vf::vec(                                                            \
            vf::exp( KOVASZNAY_LAMBDA(M_mu)*vf::Px() )*                   \
            ( (vf::pow(KOVASZNAY_LAMBDA(M_mu),2) - 4*M_PI*M_PI )          \
              *M_mu*vf::cos(2*M_PI*vf::Py())-KOVASZNAY_LAMBDA(M_mu)*vf::exp(KOVASZNAY_LAMBDA(M_mu)*vf::Px() ) ) \
            ,                                                           \
            vf::exp( KOVASZNAY_LAMBDA(M_mu)*vf::Px() )*M_mu*vf::sin(2*M_PI*vf::Py())* \
            (-vf::pow(KOVASZNAY_LAMBDA(M_mu),2) + 4*M_PI*M_PI)            \
                                                                        ) \
    /**/
#else
#define FORCE_VOLUMIQUE(NU)                                             \
    -KOVASZNAY_DIV_SIGMA(NU)                                            \
    /**/
#endif



#define CL_DEF_NS(TIME)                                                 \
    addCL<cl::dirichlet_vec,cl::ParoiB>( KOVASZNAY_VELOCITY(M_mu) )     \
    +                                                                   \
    addCL<cl::dirichlet_vec,cl::ParoiH>( KOVASZNAY_VELOCITY(M_mu) )     \
    +                                                                   \
    addCL<cl::dirichlet_vec,cl::Entree>( KOVASZNAY_VELOCITY(M_mu) )     \
    +                                                                   \
    /*addCL<cl::dirichlet_vec,cl::Sortie>( KOVASZNAY_VELOCITY(M_mu) )*/ \
    addCL<cl::neumann_scal,cl::Sortie>( KOVASZNAY_SIGMA(M_mu) )       \
    /**/

#define NS_USE_LAGMULT 0                        \
    /**/
